Force fields in GROMACS
GROMACS supports a number of force fields. From '/usr/local/gromacs/share/gromacs/top': 1: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003) 2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995) 3: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996) 4: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000) 5: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006) 6: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010) 7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002) 8: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins) 9: GROMOS96 43a1 force field 10: GROMOS96 43a2 force field (improved alkane dihedrals) 11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205) 12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656) 13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656) 14: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9) 15: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals) The Assisted Model Building with Energy Refinement (AMBER), Chemistry at Harvard Macromolecular Mechanics (CHARMM), Groningen Molecular Simulation (GROMOS), and Optimised Potentials for Liquid Simulations (OPLS) are the major families of biomolecular simulation. GROMOS fields are united-atom; the others are largely all-atom. Each family comprises a number of variants, which often perform more differently between each other than between each family. There's no standard force field for simulation, by far. A number of comparative analyses build the general conclusion that each force field is suited for particular systems, for particular tasks. These comparisons tend to be undertaken with the same simulation details, using one thermostat and barostat across the board, the same treatment of short- and long-range interactions, and the same constraints. Of course, this is not how the fields were developed -- but it is generally advisable to use the parameters that the force fields were validated against. The GROMACS site apparently used to have pages on how to deal with this using MDP options, but those have been lost to time (although CHARMM's is still up). Below are some that I used. Note that AMBER and CHARMM were developed using the Nose-Hoover thermostat, although this has been found to be nonergodic. The Nose-Hoover chain method (less nonergodic) is implemented only for the velocity Verlet (integrator = md-vv) integrator in GROMACS, which is incidentally better suited to and more accurate for pressure control simulations. Unfortunately, it requires twice as many communication calls as leapfrog (integrator = md), which slows down parallelisation. Benchmark performance on large systems before you commit to it, maybe. GROMOS Josh has a sample MDP file here. AMBER AMBER force fields tend to vary wildly in simulation parameters used. One thing to note is that the SHAKE algorithm is usually used to constrain bonds. Given that LINCS is 3-4x faster, more stable, and the same accuracy, I have chosen to stay with LINCS. Otherwise, find the paper that the force field variant was developed from for specific details. There are some considerations. Despite being an AMBER force field, the majority of papers I found using ff99sb-ildn implemented it in GROMACS. These used thermostats ranging from Nose-Hoover (the original), the Bussi-Parrinello (v-rescale in GROMACS), Langevin, and Berendsen. Barostats were usually Parrinello-Rahman, but some used Berendsen and a couple used Monte-Carlo (the only one available in AMBER). Long-range electrostatics were largely treated with PME with varying grid spacing, although some used Gaussian split Ewald. This MDP file is for ff99sb-ildn and should be applicable to most biomolecular systems for most jobs, but you should probably find a paper that's similar to what you want to do and copy them if they have good results. Note that ff99sb-ildn-nmr, ff99sb-ildn-phi, ff99sb*-ildn etc should be treated as separate force fields. Use TIP3P water. ; VARIOUS PREPROCESSING OPTIONS title = AMBER 99SB-ILDN options ; RUN CONTROL PARAMETERS integrator = md tinit = 0 ; initial time dt = 0.002 ; 2 fs nsteps = 100000000 ; 200 ns nstcomm = 1 comm-grps = Protein_Membrane Solvent ; ENERGY MINIMIZATION OPTIONS = ; Force tolerance and initial step-size = emtol = 100 emstep = 0.01 ; OUTPUT CONTROL OPTIONS nstxout = 0 ;for no .trr file nstvout = 5000 nstfout = 0 nstlog = 5000 nstenergy = 250 nstxout-compressed = 10000 compressed-x-precision = 10000 compressed-x-grps = ; groups written to xtc file energygrps = Protein Membrane Solvent ; NEIGHBORSEARCHING PARAMETERS cutoff-scheme = Verlet nstlist = 5 ; nblist update frequency, 5 x 0.002 = 10 fs ns-type = grid pbc = xyz verlet-buffer-tolerance = -1 ; override Verlet buffer tolerance for manual rlist rlist = 1.05 ; ; OPTIONS FOR ELECTROSTATICS AND VDW ; Long-range electrostatics coulombtype = PME rcoulomb-switch = 0 rcoulomb = 1.0 ; electrostatic cutoff ; Dielectric constant (DC) for cut-off or DC of reaction field = epsilon-r = 1 ; Method for doing Van der Waals vdw-type = cut-off ; vdw-modifier = Potential-shift-Verlet; rvdw-switch = 0 ; where to start switching LJ (use for force/potential switching) rvdw = 1.0 ; LJ cutoff, nm DispCorr = EnerPres ;long-range dispersion corrections fourierspacing = 0.1 fourier-nx = 0 fourier-ny = 0 fourier-nz = 0 pme-order = 4 ewald-rtol = 1e-05 epsilon-surface = 0 ; Temperature coupling tcoupl = vrescale tc-grps = Protein Membrane Solvent tau-t = 0.1 0.1 0.1 ; time constant, ps ref-t = 300 300 300 ;reference temperature, K ; Pressure coupling = Pcoupl = Parrinello-Rahman Pcoupltype = semiisotropic tau-p = 5 ; Time constant (ps) compressibility = 4.5e-5 4.5e-5 ; compressibility (1/bar) ref-p = 1.0 1.0 ; reference P (bar) ; GENERATE VELOCITIES FOR STARTUP RUN gen-vel = yes gen-temp = 300 gen-seed = -1 ;random seed ; OPTIONS FOR BONDS ; Traditionally AMBER force fields use SHAKE, but LINCS is more stable, ; faster, and as accurate. In GROMACS LINCS is normally used. constraints = h-bonds ; all bonds involving H constraint-algorithm = Lincs ; unconstrained-start = no shake-tol = 0.0001 ; for SHAKE lincs-order = 4 lincs-warnangle = 30 morse = no ; Convert harmonic bonds to morse potentials CHARMM This MDP file is for CHARMM36. I haven't investigated the other CHARMM force fields hugely, but I think that the general changes (only H-bonds constrained, force-switching, thermostat/barostat, etc) hold. Use TIP3P water. Equilibration in NPT: title = CHARMM36 eq integrator = md dt = 0.002 nsteps = 1000000 ; 2 ns nstcomm = 100 comm-grps = nstxout = 0 nstvout = 0 nstfout = 50000 nstlog = 10000 nstenergy = 10000 nstxout-compressed = 10000 compressed-x-precision = 10000 compressed-x-grps = energygrps = Protein non-Protein cutoff-scheme = Verlet nstlist = 20 ns_type = grid pbc = xyz rlist = 1.2 coulombtype = PME rcoulomb = 1.2 vdw-modifier = Force-switch rvdw_switch = 1.0 rvdw = 1.2 dispcorr = no tcoupl = Berendsen tc-grps = Protein Membrane Solvent tau_t = 1.0 1.0 1.0 ref_t = 300 300 300 Pcoupl = Berendsen Pcoupltype = semiisotropic tau_p = 1.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422 compressibility = 4.5e-5 4.5e-5 ref_p = 1.0 1.0 refcoord_scaling = COM gen_temp = 300 gen_seed = -1 constraints = h-bonds constraint_algorithm = LINCS lincs-order = 4 lincs-iter = 1 lincs-warnangle = 30 continuation = no gen_vel = yes define = -DPOSRES -DPOSRES_FC=1000 Simulation in NPT: title = CHARMM36 md integrator = md dt = 0.002 nsteps = 25000000 nstlog = 10000 nstxout = 0 nstvout = 0 nstfout = 5000 nstcalcenergy = 100 nstenergy = 10000 compressed-x-precision = 1000 nstxtcout = 10000 ; cutoff-scheme = Verlet nstlist = 20 rlist = 1.2 coulombtype = pme rcoulomb = 1.2 vdwtype = Cut-off vdw-modifier = Force-switch rvdw_switch = 1.0 rvdw = 1.2 ; tcoupl = Nose-Hoover tc_grps = Protein Membrane Solvent tau_t = 1.0 1.0 1.0 ref_t = 300 300 300 ; pcoupl = Parrinello-Rahman pcoupltype = semiisotropic tau_p = 5.0 compressibility = 4.5e-5 4.5e-5 ref_p = 1.0 1.0 ; constraints = h-bonds constraint_algorithm = LINCS continuation = no ; nstcomm = 1 comm_mode = linear comm_grps = Protein_Membrane Solvent ; refcoord_scaling = com OPLS The OPLS-UA field is rarely used for biomolecules; most common is the OPLS-AA/L one offered in GROMACS. Interestingly, OPLS was developed for Monte Carlo simulations of organic molecules in the liquid phase rather than solid phase MD. It is a bit difficult to equate conditions correctly; I followed GROMOS parameters but increased neighbour searching to 10 fs and used the Berendsen thermostat and barostat. This MDP file is for OPLS-AA/L. Dispersion corrections MUST be turned on. Use TIP4P water. title = OPLS-AA/L ; Run parameters integrator = md ; leap-frog integrator nsteps = 100000000 ; 2 * 10000000 = 20000 ps (200 ns) dt = 0.002 ; 2 fs nstcomm = 1 ; Output control ;nstxout = 0 ; save coordinates every 20 ps nstvout = 0 ; save velocities every 20 ps nstxtcout = 10000 ; xtc compressed trajectory output every 20 ps nstenergy = 10000 ; save energies every 20 ps nstlog = 10000 ; update log file every 20 ps ;nstxout-compressed = 10000 compressed-x-precision = 1000 cutoff-scheme = Verlet continuation = no ; Restarting after NPT constraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cels nstlist = 5 ; 10 fs rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm) rvdw = 1.4 ; short-range van der Waals cutoff (in nm) vdw-modifier= Potential-shift-Verlet ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.16 ; grid spacing for FFT ; Temperature coupling is on tcoupl = Berendsen ; OPLS tc-grps = Protein Membrane Solvent ; three coupling groups - more accurate tau_t = 0.1 0.1 0.1 ; time constant, in ps ref_t = 300 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Berendsen ; opls pcoupltype = semiisotropic ; uniform scaling of x-y box vectors, independent z tau_p = 0.1 ; OPLS tau ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar) compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr = EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; Velocity generation is off gen_temp = 300 gen_seed = -1 ; COM motion removal ; These options remove motion of the protein/bilayer relative to the solvent/ions comm-mode = Linear comm-grps = Protein_Membrane Solvent